Kinetic Heterogeneities in a Highly Supercooled Liquid 
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We study a highly supercooled two-dimensional fluid mixture via molecular dynamics simulation. 
We follow bond breakage events among particle pairs, which occur on the scale of the a relaxation 
time T a ■ Large scale heterogeneities analogous to the critical fluctuations in Ising systems are found 
in the spatial distribution of bonds which are broken in a time interval with a width of order 0.05r a . 
The structure factor of the broken bond density is well approximated by the Ornstein-Zernike form. 
The correlation length is of order IOOcti at the lowest temperature studied, <ti being the particle 
size. The weakly bonded regions thus identified evolve in time with strong spatial correlations. 
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As fluids are supercooled below their melting temper- 
ature, their structural relaxations become very slow, re- 
sulting in a glass transition at T = T„ characterized by 
a dramatic increase of the viscosity. |l|,|| Recently much 
attention has been paid to the mode coupling theory, 
HQ which is a self-consistent decoupling scheme for the 
density time correlation function and describes onset of 
glassy slowing down or slow structural relaxations consid- 
erably above T g . The dominant fluctuation contribution 
arises from the first peak wave number of the structure 
factor. For a long time, however, it has been intuitively 
expected (|-|^] that dynamics in glasses should be cooper- 
ative, involving many molecules, owing to configuration 
restrictions. Adam and Gibbs || speculated that particle 
motions over the interparticle distance or the potential 
barrier can take place only collectively in cooperatively 
rearranging regions (CRR) and such regions have a min- 
imum size which grows as the temperature is lowered. 
On the other hand, the mode coupling theory does not 
predict long range correlations. 

In his recent experiments Fischer || has reported 
large excess light scattering with a correlation length £ 
(20 — 200 nm) which increases on approaching the glass 
transition from a liquid state. This indicates the presence 
of large scale heterogeneities in supercooled liquids. To 
examine this effect, some authors used Monte Carlo simu- 
lations on a dense polymer melt. 1 10 llj] There, significant 
dependence of the diffusion constant on the system size 
was found, indicating correlations up to the system size. 
pot Very recently, Weber et al [jll| have found that short 
range nematic orientational order produces enhancement 
of long range density correlations. The work of Weber et 
al suggests that large scale density heterogeneities will 
not be produced in simple molecular systems. 

In their MD simulation on a simple two-dimensional 
fluid in a deeply supercooled state, Muranaka and 
Hiwatari [Q visualized significant large scale hetero- 
geneities in particle displacements in a relatively short 
time interval which corresponds to the so-called /3 relax- 
ation. In liquid states with higher temperatures, Hur- 
ley and Harrowell ]R| observed similar kinetic hetero- 
geneities but the correlation length was still on the order 



of a few particle diameters. The characterization of these 
patterns has not been made in these papers. 

It is a natural expectation that structural changes of 
the configurations of particles occur collectively M, but 
it is rather surprising that the visualization of such long 
range correlations has been nonexistent in the structural 
or a relaxation. The aims here are to visualize them 
unambiguously and analyze the patterns qualitatively. 

Our system is composed of two different atomic species, 
1 and 2, with the numbers N% — N2 = 5000, interacting 
via the soft-core potential v a /3(r) = e[(a a + ap)/2r} 12 , 
where r is the distance between two particles and a, (3 = 
1,2. We take the size and mass ratios as 02/01 = 1.4 
and rnijm\ = 2. The difference between the particle 
sizes prevents crystallization. There is no enhancement 
of the composition fluctuations at any wave number and 
the system is effectively one-component. The interac- 
tion is truncated at r = 4.50i. The leapfrog algorithm 
JLH is adopted to solve the differential equations with 
a time step of 0.005t under the periodic boundary con- 
dition. The space and time are measured in units of 
CTi and To = (midj/e) 1 / 2 . The temperature is kept at a 
constant value using the Gaussian constraint thermostat, 
|]l4| and the density is fixed at n = O.8/0 2 . We specify 
the thermodynamic states using the coupling parame- 
tcr r eJJ = n{e/k B T) 1 l 6 Y Ja ^x a xp{a a + o-pf/4, ||Jl5| 
where x\ and xi are the compositions and are equal to 
1/2 in our case, so it follows that o~ e ff — I.2I01 Data 
are taken at F e // = 1.0, 1.1, 1.2, 1.3 and 1.4. The cor- 
responding temperature is 2.54, 1.43, 0.850, 0.526 and 
0.337, respectively, in units of e/fes- The corresponding 
pressure is 32.8, 27.0, 23.6, 21.63 and 20.5, respectively, 
in units of e/0 2 . We realized these states by cooling the 
system very slowly in a stepwise manner from a liquid 
state with Y e ff — 0.8. The equilibration time is 5000 
after T e ff is changed from 1.3 to 1.4. There was no ap- 
preciable change in the pressure over computation times 
of order 10000. The zero-shear viscosity 77 is of order 6 at 
T eff = 1 and 10 4 at T ef} = 1.4 in units of mi/r @|. It 
was claimed that glassy behavior becomes apparent for 

r e// >i.3. @ 

In our analysis we follow bond breakage processes. For 
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each atomic configuration given at time to, a pair of 
atoms i and j is considered to be bonded if 



rij(t ) = \Ti(t ) - Tj(t )\ < l.la a p, 



(1) 



where i and j belong to the species a and f3, respectively, 
and <j a f3 = \(o a +<yp). In our case the particle density is 
so high that the pair distribution functions g a f}{r) have 
a sharp peak at r = o~ a p. Thus the definition of bonds 
is very insensitive to the factor in front of a a p as long 
as it is slightly larger than 1. The polygons, which are 
composed of the bonds and cover the space, are mostly 
triangles as shown in Fig. 1. After a lapse of time t 
we count the number N suv (t) of surviving bonds which 
continue to satisfy 



rij(t + t) < 1.6<7 Q(3 . 



(2) 



If the initially bonded pair does not satisfy (2), we regard 
it to have been broken. This definition of bond breakage 
is also insensitive to the factor in front of a a p as long as 
it is considerably larger than 1 and smaller than 2. We 
have confirmed that N suv (t) decays exponentially as 



AW*) = N suv (0) exp {-t/n) . 



(3) 



This holds excellently in the whole time region for T e f t < 
1.2, but has been obtained from the fit in the short time 
region t < Tfc for T e ff > 1.3. We have not yet checked 
whether or not N suv (t) relaxes exponentially for t > 77,. 
Note that there are three kinds of bonds in our two- 
component fluid, but no significant differences can be 
found in their breakage processes. In Fig. 2 we show 
that the bond breakage time tj, grows strongly with in- 
creasing T e ff. Our data may be well fitted to 77, cx T~ 4 
as well as to log 77, cx 1/T. Because 77, is the relaxation 
time of the glassy structures, 77, should be of the same 
order as the a relaxation time r a . In fact for the same 
model Muranaka and Hiwatari jlj] obtained r Q on the 
same order as our 77, from the self part of the time cor- 
relation function F s (q,t) at q = 27t/<ti. This coincidence 
is natural because F s (q,t) at q = 2w/a± decays when a 
tagged particle moves over the interparticle distance. 

Fig. 3 displays the bonds broken within the time in- 
terval [to, to + 0.0577,] at r e // = 1.0 and 1.3. The cen- 
ter positions Ry = 4(rj(fo) + Tj(to)) at the initial time 
to of the pairs broken within the time interval are de- 
picted here. For the glassy case T e ff — 1.3, the bro- 
ken bond distribution is markedly heterogeneous, being 
composed of clusters with various sizes, whereas for the 
liquid case T e ff — 1, the inhomogeneity is much weaker. 
More precisely, in the bond breakage we observe strings 
of broken bonds involving several particles even in liquid 
states, because of the high density of our system, and 
their large-scale aggregation in glassy states. The origin 
of the heterogeneities is ascribed to the fact that parti- 
cle motions over the interparticle distance can occur only 
collectively in glassy states. || Furthermore, in Fig. 4 we 
show the broken bonds in two consecutive time intervals, 



[t , t + 0.05Tb] and [t +0.05n,to + 0.1r 6 ], at T eff = 1.4. 
The initial times at which the bonds are defined are to 
and to + 0.0577, for the two cases. We can see that the 
clusters of broken bonds in the two time intervals mostly 
overlap or are adjacent to each other. This means that 
weakly bonded regions or relatively active regions migrate 
in space on the time scale of r a . The a relaxation ends 
when most of the bonds are broken. 

We next calculated the structure factor Sb(q) of the 
broken bond density pb(v) = ^ S(r — R, 3 ) defined by 



Sb(q) = N b 



y%xp(»q- Rjj) 



(4) 



The summation is over the broken pairs, Nb is the to- 
tal number of the broken bonds, and (•••)„ is the an- 
gular average over the direction of the wave vector q = 
(2ir/ L)(n x , n y ) where n x , n y — ±1, ±2, ±3. • • •. The sys- 
tem length L is 111.8 here. Furthermore, since Sb (q) from 
one configuration fluctuates at small q, we have taken 
the averages of Sb(q) from 10 sequential configurations 
for T e ff = 1.0 — 1.3 and 4 sequential configurations for 
T e ff — 1-4. Then, as shown in Fig. 5, Sb{q) can be well 
approximated by the simple Ornstcin-Zernike form, 



S b (q)=S b (0)/[l + eq% 



(5) 



from which we can determine £. In addition, the raw 
Sb(q) from (4) tends to become 1 at large q, so it has 
been subtracted in Sb(q) in Fig. 5. (This subtracted 
value is negligibly small for T e ff > 1.1 at long wave- 
lengths qt; < 1, however.) Interestingly, 



s b (o) cx e 



(6) 



and the large q behavior of Sb {q) is insensitive to T e f f as 
in Ising spin systems near the critical point. In Fig. 6 we 
show £ vs 1/T for the two widths of the time intervals, 
0.05rfc (O) and O.ln (•). We can see that £ is largely 
independent of the width except for the lowest temper- 
ature case, where T e ff = 1-4. In this case, however, £ 
approaches the system length L — 111.8 and a finite size 
effect should be present. Also 77, at the lowest T in Fig. 
2 seems to be affected by the same effect. 

Figures 3-5 suggest that the cluster structure in the 
highly supercooled states comprises clusters with various 
length scales, the minimum being o~\ and the maximum 
being £. It appears to be self-similar if we consider clus- 
ters with size I in the region o\ <C t <§C £ in the limit 
£ 3> <ti, which is also supported by the power law be- 
havior Sb(q) ~ 1/<Z 2 shown in Fig. 5. This apparently 
contradicts Adam-Gibbs' phenomenological speculation 
in which the minimum size of CRR grows upon cooling 
||. We stress that the heterogeneities are of purely ki- 
netic origin and are dynamical objects. Indeed, the small 
clusters with I <C £ in Figs. 3 and 4 evolve into larger 
ones or dissolve in the subsequent time intervals. 

We have also calculated the structure factors of the 
density, composition and stress tensor, but no indica- 
tion of large scale heterogeneities has been detected in 
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such static quantities in accordance with previous pa- 
pers, (ljjl^] The kinetic heterogeneities found in this pa- 
per play a crucial role in response to externally applied 
perturbations such as sound waves or shear flow. ||16| 

A frequently disputed issue in the literature is whether 
or not there is an underlying thermodynamic phase tran- 
sition at a nonzero temperature To below T g in highly 
supercooled liquids. In our case, Fig. 6 does not indicate 
any divergence of £ at a nonzero temperature, although 
this is not conclusive due to the finite size effect aris- 
ing from £ ~ L. We mention a recent simulation on a 
three-dimensional Lennard- Jones mixture by Ghosh and 
Dasgupta, |2(J who used scaling relations of finite sys- 
tems to conclude £ocT~ 15 asT— > without examining 
cluster structures. The limited data in Fig. 6 also yield 
|<91og£/<91ogT| ~ 1. 

In summary, by performing long time and large scale 
MD simulations on a 2D soft-core mixture, we have ex- 
amined bond breakage events to find large-scale hetero- 
geneities in the changing rates of the amorphous struc- 
ture, which are enhanced at low temperatures and sur- 
prisingly similar to the critical fluctuations in Ising sys- 
tems. We believe that this is a first step in the inves- 
tigation of more complex situations and 3D cases. We 
should also clarify the relationship of our patterns in the 
a relaxation and those by Muranaka and Hiwatari fl2| | 
in the (3 relaxation. 
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FIGURE CAPTIONS 

Fig.l Bonds defined at a given time. Their lengths are 
mostly close to a a p. 

Fig. 2 Temperature dependence of the bond breakage 
time Tb (•), where T* — ksT/e. The dotted line is a 
viewing guide. The vertical arrows are estimates for the 
a relaxation time from the MD data of Muranaka and 
Hiwatari. Il7j 

Fig. 3 Snapshots of the broken bonds at T e ff = 1.0 and 
1.3. The time interval is 0.05rfc, so 5% of the initial bonds 
are broken here. The arrows indicate £. 

Fig. 4 Broken bond distributions in two consecutive time 
intervals, [t ,t + 0.05r 6 ] (□) and [t + 0.05r fc , t + O.ln] 
(•), at T e ff = 1.4. The arrows indicate £. 

Fig. 5 Sb(q) of the broken bond density. The insert shows 
1/Sb(q) vs g 2 , from which £ -2 is determined. 

Fig. 6 Temperature dependence of the correlation length 
£, where T* — ksT/e. The widths of the time intervals 
are 0.05t{, (O) and O.lr;, (•). The system length is 111.8. 
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